MODULE measure
  !-----------------------------------------------------------
  ! Measures of goodness used by the KRUUNUNHAKA basis set tool kit
  !-----------------------------------------------------------
  ! See CHANGELOG
  !-----------------------------------------------------------
  ! Copyright (C) 2003-2008 Pekka Manninen, 2010 Jussi Lehtola
  !
  ! This program is distributed under the terms of the GNU General
  ! Public License as published by the Free Software Foundation;
  ! either version 2 of the License, or (at your option) any later version.
  !
  ! This program is distributed in the hope that it will be useful,
  ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
  !
  ! Please note that this program as well as its 
  ! author must be properly cited whenever the program or some parts 
  ! originated on it are employed. 
  !
  !-----------------------------------------------------------
  
  USE definitions
  USE completeness_profile

  IMPLICIT NONE
CONTAINS

  SUBROUTINE f_everything(l, X, C, alphas, usesto)
    ! Calculate f_dev, f_lsq and f_sum all at once and print out result

    ! Angular momentum value
    INTEGER, INTENT(IN) :: l
    ! Exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:), INTENT(IN) :: C
    ! Values of exponentials to scan with
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Scan with STOs?
    LOGICAL, INTENT(IN) :: usesto

    ! Measures of goodness
    REAL(KIND=prec) :: maxdev, lsq, sum
    ! Logarithms of scanning exponents
    REAL(KIND=prec), DIMENSION(SIZE(alphas)) :: loga
    ! Completeness profile
    REAL(KIND=prec), DIMENSION(SIZE(alphas)) :: Y

    ! Helpers
    REAL(KIND=prec) :: ldiff, mdiff, rdiff
    ! Loop index
    INTEGER :: i

    ! Compute log(alpha)
    loga=LOG(alphas)
    ! Compute profile
    CALL comp_prof(l, X, C, alphas, usesto, Y)

    ! Maximum deviation from unity
    maxdev=MAXVAL(1.0_prec-Y)

    ! Compute integrals
    lsq=0.0_prec
    sum=0.0_prec
    DO i=2, SIZE(alphas)-1, 2
       ! Compute left, middle and right value
       ldiff = 1.0_prec - Y(i-1)
       mdiff = 1.0_prec - Y(i)
       rdiff = 1.0_prec - Y(i+1)
       ! Increment lsq and sum
       lsq = lsq + (ldiff**2 + 4.0_prec*mdiff**2 + rdiff**2)/6.0_prec*(loga(i+1)-loga(i-1))
       sum = sum + (ldiff + 4.0_prec*mdiff + rdiff)/6.0_prec*(loga(i+1)-loga(i-1))
    END DO
    ! Normalize wrt length of interval
    lsq = SQRT(lsq / (loga(SIZE(alphas))-loga(1)))
    sum = sum / (loga(SIZE(alphas))-loga(1))

    WRITE(*,'(A50,G11.4,/,A50,G11.4,/,A50,G11.4,/A50,/)') &
         'Maximum deviance from completeness:', maxdev, &
         'Integrated least-squares difference:', lsq, &
         'Integration over the profile:', 1.0_prec-sum, &
         '(All normalized to the length of the interval)'

  END SUBROUTINE f_everything

    
  FUNCTION f_maxdev(l, X, C, alphas, usesto)
    ! Maximum deviance from unity

    ! Angular momentum value
    INTEGER, INTENT(IN) :: l
    ! Exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:), INTENT(IN) :: C
    ! Values of exponentials to scan with
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Scan with STOs?
    LOGICAL, INTENT(IN) :: usesto

    ! Completeness profile
    REAL(KIND=prec), DIMENSION(SIZE(alphas)) :: Y

    ! Result
    REAL(KIND=prec) :: f_maxdev

    ! Compute profile
    CALL comp_prof(l, X, C, alphas, usesto, Y)

    f_maxdev = MAXVAL(1.0_prec - Y)
  END FUNCTION f_maxdev


  FUNCTION f_sum(l, X, C, alphas, usesto)
    ! Summed deviation from unity / integral over profile

    ! Angular momentum value
    INTEGER, INTENT(IN) :: l
    ! Exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:), INTENT(IN) :: C
    ! Values of exponentials to scan with
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Scan with STOs?
    LOGICAL, INTENT(IN) :: usesto

    ! Logarithms of scanning exponents
    REAL(KIND=prec), DIMENSION(SIZE(alphas)) :: loga
    ! Completeness profile
    REAL(KIND=prec), DIMENSION(SIZE(alphas)) :: Y

    ! Loop index
    INTEGER :: i

    ! Result
    REAL(KIND=prec) :: f_sum

    ! Compute log(alpha)
    loga=LOG(alphas)
    ! Compute profile
    CALL comp_prof(l, X, C, alphas, usesto, Y)

    ! Compute integral over profile
    f_sum=0.0_prec
    DO i=2, size(alphas)-1, 2
       f_sum = f_sum + 0.5_prec*(Y(i-1)+4.0_prec*Y(i)+Y(i+1))*(loga(i+1)-loga(i-1))
    END DO
    ! Normalize with respect to length of integration interval
    f_sum = f_sum / (loga(SIZE(alphas))-loga(1))
    ! and finally convert the result into a deviation from completeness
    f_sum = 1.0_prec - f_sum
  END FUNCTION f_sum


  FUNCTION f_lsq(l, X, C, alphas, usesto)
    ! Sum of squares of differences from unity

    ! Angular momentum value
    INTEGER, INTENT(IN) :: l
    ! Exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:), INTENT(IN) :: C
    ! Values of exponentials to scan with
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Scan with STOs?
    LOGICAL, INTENT(IN) :: usesto
    
    ! Logarithms of scanning exponents
    REAL(KIND=prec), DIMENSION(SIZE(alphas)) :: loga
    ! Completeness profile
    REAL(KIND=prec), DIMENSION(SIZE(alphas)) :: Y

    ! Loop index
    INTEGER :: i

    ! Helpers
    REAL(KIND=prec) :: ldiff, mdiff, rdiff

    ! Result
    REAL(KIND=prec) :: f_lsq

    ! Compute log(alpha)
    loga=LOG(alphas)
    ! Compute profile
    CALL comp_prof(l, X, C, alphas, usesto, Y)

    f_lsq=0.0_prec
    ! Compute integral using 3-point Simpson rule
    DO i=2, SIZE(alphas)-1, 2
          ! Compute left, middle and right value
          ldiff = (1.0_prec - Y(i-1))
          mdiff = (1.0_prec - Y(i))
          rdiff = (1.0_prec - Y(i+1))
          ! Increment f_lsq
          f_lsq = f_lsq + (ldiff**2 + 4.0_prec*mdiff**2 + rdiff**2)/6.0_prec*(loga(i+1)-loga(i-1))
    END DO

    ! Normalize wrt length of interval
    f_lsq = SQRT(f_lsq / (loga(SIZE(alphas))-loga(1)))
  END FUNCTION f_lsq
END MODULE measure

  
    

  
